The response to selection, sexual antagonism and sexual concordance in the DGRP

workflowr

• Summary
• Checks
• Past versions

Last updated: 2024-04-11

Checks: 7 0

Knit directory: DGRP_sexual_conflict/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20210706)

The command set.seed(20210706) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 3181c8d

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 3181c8d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rapp.history
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  %
    Untracked:  Antagonism_trait_plot.pdf
    Untracked:  Chapter_2_Figure_S1.pdf
    Untracked:  Chapter_2_Figure_S2.pdf
    Untracked:  Conflict_plot_females.pdf
    Untracked:  Conflict_plot_male.pdf
    Untracked:  Diff_atlas_plot.pdf
    Untracked:  Figure_3.pdf
    Untracked:  Figure_4.pdf
    Untracked:  Figures.Rmd
    Untracked:  Mainz_response_space_plot.pdf
    Untracked:  Manuscript/
    Untracked:  Morrissey_model.html
    Untracked:  Morrissey_model.qmd
    Untracked:  Mutation_load_phenotype_effect.Rmd
    Untracked:  PRISMA.pptx
    Untracked:  R_atlas_plot.pdf
    Untracked:  R_multivariate_brms_simple.rds
    Untracked:  R_trans_multivariate_brms_simple.rds
    Untracked:  R_transcriptome_medians.csv
    Untracked:  Reported_heritability.xlsx
    Untracked:  Rplot.pdf
    Untracked:  Rplot01.tiff
    Untracked:  Useful_cuts.Rmd
    Untracked:  add_later.Rmd
    Untracked:  antagonism_R_space_plot.svg
    Untracked:  antagonism_trait_plot.svg
    Untracked:  big_f_mutation_plot.pdf
    Untracked:  big_m_mutation_plot.pdf
    Untracked:  code/get_gene_annotations.R
    Untracked:  conflict_plot.pdf
    Untracked:  data/DGRP_mutation_load.csv
    Untracked:  data/R_medians.csv
    Untracked:  data/X_burden.csv
    Untracked:  data/autosomal_burden.csv
    Untracked:  data/female_X_mutation_corrs.csv
    Untracked:  data/female_autosomal_mutation_corrs.csv
    Untracked:  data/female_mutation_corrs.csv
    Untracked:  data/input/
    Untracked:  data/male_X_mutation_corrs.csv
    Untracked:  data/male_autosomal_mutation_corrs.csv
    Untracked:  data/male_mutation_corrs.csv
    Untracked:  data/organismal_level_output/
    Untracked:  data/transcriptome_output/
    Untracked:  evidence_ratios.R
    Untracked:  female_mutation_load.pdf
    Untracked:  figures/
    Untracked:  fits/
    Untracked:  height_example.pdf
    Untracked:  mag_R_sem.pdf
    Untracked:  male_mutation_load.pdf
    Untracked:  mutation_load_plot.pdf
    Untracked:  sem_figures.R
    Untracked:  sim_me.rds
    Untracked:  sim_no_me.rds

Unstaged changes:
    Modified:   _workflowr.yml
    Modified:   analysis/Transcriptome_analysis.Rmd
    Deleted:    code/export_gwas_results_for_shiny_app.R
    Deleted:    data/all.dgrp.phenos_unscaled.csv

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Main_analysis.Rmd) and HTML (docs/Main_analysis.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	3181c8d	tomkeaney	2024-04-11	Completed analysis
html	24fdbb4	tomkeaney	2024-04-11	Build site.
Rmd	1228613	tomkeaney	2024-04-11	Completed analysis
html	1228613	tomkeaney	2024-04-11	Completed analysis
html	ad196cf	tomkeaney	2024-04-11	Build site.
Rmd	ec2fc6d	tomkeaney	2024-04-11	Completed analysis??
html	5b2cdee	tomkeaney	2023-10-19	Build site.
Rmd	4188afb	tomkeaney	2023-10-19	small wording changes
html	319c84c	tomkeaney	2023-10-17	Build site.
Rmd	48c7e9f	tomkeaney	2023-10-17	Restarting the project
html	92d2a26	ausevo	2023-03-08	Build site.
Rmd	5657ce7	ausevo	2023-03-08	Better Table
html	414b217	ausevo	2023-03-08	Build site.
Rmd	2501948	ausevo	2023-03-08	Change to Tom’s analysis
html	c33ab4e	ausevo	2023-03-08	Build site.
Rmd	ca66d6a	ausevo	2023-03-08	Change to Tom’s analysis
html	daa935d	ausevo	2023-03-08	Build site.
Rmd	78ece4f	ausevo	2023-03-08	Change to Tom’s analysis
html	8f6bea5	ausevo	2023-03-03	Build site.
Rmd	a382ebf	ausevo	2023-03-03	New version
html	b682a4f	ausevo	2023-03-03	Build site.
Rmd	6577470	ausevo	2023-03-03	New version
Rmd	7bb248d	ausevo	2023-01-31	rework of line means results
html	6d630f0	ausevo	2022-05-12	Build site.
Rmd	95f4a97	ausevo	2022-05-12	including ideas section
html	7bf0ee0	ausevo	2022-05-10	Build site.
Rmd	5b9883e	ausevo	2022-05-10	Add GCTA model outputs
Rmd	8e2d8be	tkeaney	2021-08-17	Merge branch ‘master’ of https://github.com/tomkeaney/DGRP_sexual_conflict
Rmd	4b095f0	tkeaney	2021-08-17	edits
Rmd	73f011c	lukeholman	2021-08-06	Tiny edit (digits=3)
html	043902f	tkeaney	2021-07-28	Build site.
Rmd	959123f	tkeaney	2021-07-28	trying to fix a bug
html	769e48d	tkeaney	2021-07-28	Build site.
Rmd	340c41b	tkeaney	2021-07-28	sexual dimorphism calculation added
html	4018580	tkeaney	2021-07-28	Build site.
Rmd	65b2597	tkeaney	2021-07-28	sexual dimorphism calculation added
html	83a0c4f	tkeaney	2021-07-28	Build site.
Rmd	ddea55e	tkeaney	2021-07-28	sexual dimorphism calculation added
html	b364d30	tkeaney	2021-07-23	Build site.
Rmd	e404cf9	tkeaney	2021-07-23	progression on selection calculations
html	32f107f	tkeaney	2021-07-23	Build site.
Rmd	6aad8a8	tkeaney	2021-07-23	progression on selection calculations
html	11a8390	tkeaney	2021-07-07	Build site.
Rmd	42b8f12	tkeaney	2021-07-07	Get the site up and running

---

Load packages and the data

\(~\)

First load the packages and build helper functions

library(tidyverse) # for tidy coding
library(MetBrewer) # for many nice colour palettes
library(rcartocolor) # more cool colours
library(kableExtra) # for scrolling tables
library(DT) # for interactive tables
library(patchwork) # to join multiple plots nicely
library(brms) # for bayesian models
library(tidybayes) # for more bayesian things
library(bayestestR) # for the pd metric 
library(broom) # convert results of functions into tables
library(ggtext) # for markdown features in ggplot
library(ggrepel) # for plot labels in ggplot
library(ggnewscale) # to reset scales in ggplot 
library(bigsnpr) 

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      autoWidth = TRUE,
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=1000,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'Trait_data')),
      pageLength = 100
    )
  )
}

\(~\)

Organismal-level trait analysis

Selecting data appropriate for analysis

We conducted a near-exhaustive search of the literature until January 2022, to obtain line mean estimates and associated meta-data for quantitative traits that have been measured in the DGRP. We did not include data collected for traits that had been measured in heterozygous combinations of multiple DGRP lines. In total, we identified 125 studies that reported line means or raw data for 2115 phenotypic traits. To ready the data for analysis, we grouped trait values by trait and sex and standardised the data to have mean = 0 and sd = 1. Sex-specific standardised line means for each trait were then combined with their standardised fitness estimates, obtained from Wong and Holman (2023). We also include some helpful metadata for downstream analysis. We then pruned the dataset to only include traits that have been measured in single-sex cohorts, in 80 or more lines. We also removed traits included in two large datasets on the microbiome and metabolome (Everett et al. 2020 and Jin et al. 2020).

# load in the data, note that traits have already been standardised

DGRP_data <- 
  left_join(
    read_csv("data/input/all.dgrp.phenos_scaled.csv") %>% 
      mutate(line = as.factor(line)),
    read_csv("data/input/meta_data_for_all_traits.csv") %>%
      group_by(Reference) %>% 
      mutate(study_ID = as.factor(cur_group_id()),
             Pooled = if_else(Sex == "Pooled", "Yes", "No"))) %>% 
  left_join(read_rds("data/input/trait_names.rds"))
  
  
# Apply the selection criteria with the filter() function, then add a column each for female and male fitness

trait_data <-
  DGRP_data %>%  
  filter(!str_detect(Trait, "fitness"),
         `# lines measured` >= 80 &
           Pooled != "Yes" & 
           Reference != "Jin et al (2020) PLOS Genetics" & 
           Reference != "Everett et al (2020) Genome Research" &
           !str_detect(Trait, "dopamine.response.to.paraquat.2021.m")) %>% # data for this trait was entered into the database incorrectly, it was only measured in 10 lines so should not be included in our analysis 
  
  # join the early life fitness data from Wong and Holman
  
  left_join(
    DGRP_data %>%
      filter(str_detect(Trait, "fitness.early.life")) %>% 
      select(line, Trait, trait_value) %>% 
      pivot_wider(names_from = Trait, values_from = trait_value) %>% 
      rename(female_fitness = fitness.early.life.f, male_fitness = fitness.early.life.m))


clean_meta_data <-
  trait_data %>% 
  select(Trait_nice, Trait, Life_stage, `Trait guild`, study_ID, Trait_nice, Reference, `Trait description`) %>% 
  distinct(Trait, .keep_all = TRUE) %>% 
    mutate(Trait_nice = case_when(Trait_nice == "DDT resistance mortality" ~ "DDT susceptibility",
                                  .default = Trait_nice)) # make confusing name simpler

\(~\)

Table S1. Traits included for analysis

# how many studies did we start with?

table_data <- left_join(
  DGRP_data %>%
    distinct(Trait, .keep_all = TRUE) %>% 
    mutate(Measured_in = case_when(Sex == "Female" ~ "Females",
                                   Sex == "Male" ~ "Males",
                                   Sex == "Pooled" ~ "Mixed sex cohorts")) %>% 
    select(Trait_nice, Trait, Reference, `# lines measured`, Measured_in, `Trait description`),
  
  trait_data %>%
    distinct(Trait) %>% 
    mutate(Included = "Yes")
) %>%
  mutate(Included = if_else(is.na(Included), "No", Included)) %>% 
  select(-Trait) %>% 
  rename(Trait = Trait_nice, `Measured in` = Measured_in) %>% 
  arrange(desc(Included))

# Create a function to build HTML searchable tables

my_data_table(table_data %>% 
               select(Trait, Reference, `# lines measured`, `Measured in`, Included, `Trait description`))

\(~\)

Find the number of traits and studies included in our analysis.

# how many studies and traits do we have after filtering?

num_unique_traits <- table_data %>% filter(Included == "Yes") %>% nrow()
# in females
num_unique_traits_f <- table_data %>% filter(Included == "Yes" & `Measured in` == "Females") %>% nrow()
# in males
num_unique_traits_m <- table_data %>% filter(Included == "Yes" & `Measured in` == "Males") %>% nrow()
# how many studies are they measured across in total?
num_unique_studies <- table_data %>% filter(Included == "Yes") %>%  distinct(Reference) %>% nrow()
# in females
num_unique_studies_f <- table_data %>% filter(Included == "Yes" & `Measured in` == "Females") %>% distinct(Reference) %>% nrow()
# in males
num_unique_studies_m <- table_data %>% filter(Included == "Yes" & `Measured in` == "Males") %>% distinct(Reference) %>% nrow()

After this selection process, 474 remain, that were measured across 76 studies. There are 232 measured in females and 242 in males across 56 and 54 respectively.

\(~\)

Calculating R: the response to selection

\(~\)

We assume that line represents a genetically homogeneous, homozygous genotype, such that the line mean for a particular trait equals the breeding value for that trait for that genotype. We can therefore estimate the genetic variance for a given trait in the DGRP as the variance in line means (Mackay, 2012).

We can then estimate the response to selection (\(R\)) for a trait using Robertson’s Secondary Theorem of Natural Selection (also known as the Robertson covariance; Robertson, 1968), which states that \(R\) for a given trait is equivalent to the additive genetic covariance between the trait (\(A_z\)) and relative fitness (\(A_w\)), such that \(R = \sigma(A_w, Az)\). We can estimate this genetic covariance as the correlation in breeding values for traits with relative fitness across DGRP lines.

In species with two sexes and obligate sexual reproduction, half of the response to selection is due to selection on females and half to selection on males (e.g. Lande, 1980). We can therefore partition the secondary theorem into the response via each sex:

\[R = \frac{\sigma(A_w^F, A_z) + \sigma(A_w^M, A_z)}{2}\]

Furthermore, in quantitative genetics it is common to regard the same trait expressed in males or females (e.g. female body size and male body size) as separate traits. For convenience in the code below, we denote the response to selection of a female trait (\(R_F\)) or a male trait (\(R_M\)) as:

\[R_F = \frac{R_{FF} + R_{MF}}{2}\]

\[R_M = \frac{R_{MM} + R_{FM}}{2}\]

Where \(R_{FF}\) is the response of a female trait to selection on females, \(R_{MF}\) is the response of a female trait to selection on males, \(R_{MM}\) is the response of a male trait to selection on males, and \(R_{FM}\) is the response of a male trait to selection on females. These four quantities can be found via the partitioned Robertson covariance as follows:

\[R_{FF} = \sigma(A_w^F, A_z^F)\]

\[R_{MF} = \sigma(A_w^M, A_z^F)\]

\[R_{MM} = \sigma(A_w^M, A_z^M)\]

\[R_{FM} = \sigma(A_w^F, A_z^M)\]

Note that even though female traits are not expressed in males and therefore are not under selection in males, selection on males can still generate a response to selection in female traits if there is genetic covariance between these female traits and male traits that are under selection in males (and similarly, male traits can respond to selection on females due to between-sex genetic correlations). We therefore expect \(R_{MF}\) and \(R_{FM}\) to commonly be non-zero.

\(~\)

Build models to calculate \(R_{FF}\) & \(R_{MF}\)

To estimate the covariance between \(A_w\) and \(A_z\), we fitted bivariate Bayesian linear models using the brms package (Bürkner, 2017) for R version 4.2.2. For each combination of trait and sex, we used line means for the focal trait and the fitness of the focal sex as the two response variables and fitted an intercept-only Gaussian model. Each model returned a posterior distribution of the residual correlation between trait and fitness, which for data expressed in standard units is equivalent to the covariance.

Build functions to run the models

# RFF estimates

female_traits <- trait_data %>% filter(Sex == "Female")

trait_list_female <- unique(female_traits$Trait) # an input to the map_dfr() function that we'll need in a few chunks time

# code the model structure we will use for all traits using one example - `flight.performance.f`. We can then use the update() function to run this model many times, once for each trait measured in females. update() makes this process many times faster, because the model can immediately start sampling, without the need to recompile. 

RFF_model <-
  brm(data = female_traits %>% filter(Trait == "flight.performance.f"),
      family = gaussian,
      bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = femalefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1, file = "fits/RFF_test_model")   

# make a function to update the model and the posterior sample output with the 'selected trait'

RFF_calculator <- function(selected_trait){
  
  data <- female_traits %>% filter(Trait == selected_trait)
  
  model <- update(
    RFF_model, newdata = data,
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>% 
    mutate(Trait = selected_trait) %>% 
    select(Trait, Response_to_selection_female) %>% 
    as_tibble()
  
  posterior
}

# RMF estimates

RMF_model <-
  brm(data = female_traits %>% filter(Trait == "flight.performance.f"),
      family = gaussian,
      bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = malefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1, file = "fits/RMF_test_model")

# make a function to update the model and the posterior sample output with your desired trait

RMF_calculator <- function(selected_trait){
  
  data <- female_traits %>% filter(Trait == selected_trait)
  
  model <- update(
    RMF_model, newdata = data,
    chains = 4, cores = 4, iter = 6000, warmup = 2000,
    seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>% 
    mutate(Trait = selected_trait) %>% 
    select(Trait, Response_to_selection_male) %>% 
    as_tibble()
  
  posterior
}

\(~\)

Build models to calculate \(R_{FM}\) & \(R_{MF}\)

\(~\)

# RMM estimates

male_traits <- trait_data %>% filter(Sex == "Male")

trait_list_male <- unique(male_traits$Trait)

RMM_model <-
  brm(data = male_traits %>% filter(Trait == "flight.performance.m"),
      family = gaussian,
      bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = malefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1, file = "fits/RMM_test_model")   

# make a function to update the model and the posterior sample output with your desired trait

RMM_calculator <- function(selected_trait){
  
  data <- male_traits %>% filter(Trait == selected_trait)
  
  model <- update(
    RMM_model, newdata = data,
    chains = 4, cores = 4, iter = 6000, warmup = 2000,
    seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
    mutate(Trait = selected_trait) %>% 
    select(Trait, Response_to_selection_male) %>% 
    as_tibble()
  
  posterior
}

# RFM estimates

RFM_model <-
  brm(data = male_traits %>% filter(Trait == "flight.performance.m"),
      family = gaussian,
      bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = femalefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1, file = "fits/RFM_test_model")   

# make a function to update the model and the posterior sample output with your desired trait

RFM_calculator <- function(selected_trait){
  
  data <- male_traits %>% filter(Trait == selected_trait)
  
  model <- update(
    RFM_model, newdata = data,
      chains = 4, cores = 4, iter = 6000, warmup = 2000,
      seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>% 
    mutate(Trait = selected_trait) %>% 
    select(Trait, Response_to_selection_female) %>% 
    as_tibble()
  
  posterior
}

\(~\)

Run the models for all the traits

Run the models using RFF_calculator, RMF_calculator, RMM_calculator and RFM_calculator

# run the RFF function

Run_function <- FALSE # change this to TRUE to run the models

if(Run_function){
  RFF <- map_dfr(trait_list_female, RFF_calculator) # map_dfr returns a data frame created by row-binding each output
  write_csv(RFF, file = "data/organismal_level_output/RFF.csv")
} else RFF <- read_csv("data/organismal_level_output/RFF.csv")

# run the RMF function

if(Run_function){
  RMF <- map_dfr(trait_list_female, RMF_calculator) 
  write_csv(RMF, file = "data/organismal_level_output/RMF.csv")
} else RMF <- read_csv("data/organismal_level_output/RMF.csv")

# run the RMM function

if(Run_function){
  RMM <- map_dfr(trait_list_male, RMM_calculator) 
  write_csv(RMM, file = "data/organismal_level_output/RMM.csv")
} else RMM <- read_csv("data/organismal_level_output/RMM.csv")

# run the RFM function

if(Run_function){
  RFM <- map_dfr(trait_list_male, RFM_calculator) 
  write_csv(RFM, file = "data/organismal_level_output/RFM.csv")
} else RFM <- read_csv("data/organismal_level_output/RFM.csv")

\(~\)

Plot the sex-specific responses to selection

R_female_traits <-
  bind_cols(RFF, RMF) %>% 
  rename(Trait = Trait...1) %>% 
  mutate(Trait_Sex = "Female") %>% 
  select(-(Trait...3))

p1 <-
  R_female_traits %>% 
  group_by(Trait) %>% 
  mutate(avg_R = median(Response_to_selection_female)) %>%
  ggplot(aes(Response_to_selection_female, fct_reorder(Trait, avg_R))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Purp") +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
  labs(x = "Response to selection in females",
       y = "Trait expressed in females") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))

p2 <-
  R_female_traits %>% 
  group_by(Trait) %>% 
  mutate(avg_R = median(Response_to_selection_male)) %>%
  ggplot(aes(Response_to_selection_male, fct_reorder(Trait, avg_R))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Peach") +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
  labs(x = "Response to selection in males",
       y = "Trait expressed in females") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8)) 

p1 + p2 +
  plot_annotation(tag_levels = 'a')

Past versions of unnamed-chunk-7-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11
319c84c	tomkeaney	2023-10-17
daa935d	ausevo	2023-03-08
b682a4f	ausevo	2023-03-03
7bf0ee0	ausevo	2022-05-10
769e48d	tkeaney	2021-07-28
32f107f	tkeaney	2021-07-23
11a8390	tkeaney	2021-07-07

Figure S1. The estimated response to selection in a females and b males for all traits measured in females. Innermost bands approximate the median, while outer bands show the 66 and 95% credible intervals.

\(~\)

R_male_traits <-
  bind_cols(RFM, RMM) %>% 
  rename(Trait = Trait...1) %>%  
  mutate(Trait_Sex = "Male") %>% 
  select(-(Trait...3))

p3 <-
  R_male_traits %>% 
  group_by(Trait) %>% 
  mutate(avg_R = median(Response_to_selection_male)) %>%
  ggplot(aes(Response_to_selection_male, fct_reorder(Trait, avg_R))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Peach") +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
  labs(x = "Response to selection in males",
       y = "Trait expressed in males") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))

p4 <- 
  R_male_traits %>% 
  group_by(Trait) %>% 
  mutate(avg_R = median(Response_to_selection_female)) %>%
  ggplot(aes(Response_to_selection_female, fct_reorder(Trait, avg_R))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Purp") +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
  labs(x = "Response to selection in females",
       y = "Trait expressed in males") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8)) 

p3 + p4 +
  plot_annotation(tag_levels = 'a')

Past versions of unnamed-chunk-8-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11
b682a4f	ausevo	2023-03-03
7bf0ee0	ausevo	2022-05-10
043902f	tkeaney	2021-07-28
83a0c4f	tkeaney	2021-07-28
32f107f	tkeaney	2021-07-23

Figure S2. The response to selection in a males and b females for all traits measured in males. Innermost bands approximate the median, while outer bands show the 66 and 95% credible intervals.

Calculate various metrics

\(~\)

Combine the data frames and estimate the overall expected response to selection as

\[R = \frac{R_{F} + R_{M}}{2}\]

In the same code chunk, for each trait, we also calculate the difference in the response to selection acting on females and males as \(\Delta R = R_{FF} - R_{MF}\) for female traits or \(\Delta R = R_{FM} - R_{MM}\) for male traits. This difference facilitates identification of traits where the response to selection is very different between sexes (in magnitude, and possibly also in sign). Note that \(\Delta R\) alone does not reveal whether a trait is sexually antagonistic, because sexually concordant selection that varies in strength between sexes (or sex-limited selection) can result in a high value of \(\Delta R\).

R_female_traits <-
  R_female_traits %>% 
  mutate(R_overall = (Response_to_selection_female + Response_to_selection_male)/2,
         R_diff = Response_to_selection_female - Response_to_selection_male) %>% 
    select(Trait, Trait_Sex, everything()) 

R_male_traits <-
  R_male_traits %>%  
  mutate(R_overall = (Response_to_selection_female + Response_to_selection_male)/2,
         R_diff = Response_to_selection_female - Response_to_selection_male) %>% 
  select(Trait, Trait_Sex, everything()) %>% 
  filter(Trait != "dopamine.response.to.paraquat.2021.m")

R_all_traits <- bind_rows(R_female_traits, R_male_traits)

R_long_form <-
  R_all_traits %>% 
  select(1:4) %>% 
    pivot_longer(cols = 3:4, names_to = "Fitness_Sex", values_to = "R_metric") %>% 
    mutate(Fitness_Sex = case_when(Fitness_Sex == "Response_to_selection_female" ~ "Female",
                                   Fitness_Sex == "Response_to_selection_male" ~ "Male"))

\(~\)

Calculating evidence ratios

We use the brms hypothesis() function to compute evidence ratios, where our hypothesis is \(R\) (and its various related metrics) != 0. When the hypothesis is one-sided, this is the posterior probability under the hypothesis against its alternative. That is, when the hypothesis is of the form R > 0, the evidence ratio is the ratio of the posterior probability of R > 0 and the posterior probability of R < 0. In this example, values greater than one indicate that the evidence in favour of R > 0 is larger than evidence in favour of R < 0. In contrast, values smaller than one indicate that there is greater evidence in favour of R < 0 than R > 0. More on the hypothesis function can be found here. We consider evidence ratios > 39 or < 1/39 as biologically notable (without accounting for multiple testing), as our hypothesis is two-directional (negative and positive values of \(R\) are of interest to us) and these values closely correspond to the familiar frequentist p-value = 0.05 (Makowski et al, 2019).

# build a function to calculate the evidence ratio, get the relevant info from the output and convert it to a tibble

find_evidence_ratios <- function(Trait_name, specify_hypothesis){
  x <- hypothesis(R_all_traits %>% 
                    filter(Trait == Trait_name),
                  specify_hypothesis)
  
  x <- x$hypothesis
  
  x %>% as_tibble() %>% 
    select(Evid.Ratio) %>% 
    mutate(Trait = Trait_name)
}

# list of traits we need evidence ratios for.

all_traits_list <- R_all_traits %>% distinct(Trait) %>% pull(Trait)

# calculate evidence ratios 

if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){

R_evidence_ratios_female <- 
  map_dfr(all_traits_list, 
          find_evidence_ratios, 
          "Response_to_selection_female > 0") %>% 
  rename(Female_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_male <-
  map_dfr(all_traits_list, 
          find_evidence_ratios, 
          "Response_to_selection_male > 0") %>% 
  rename(Male_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_overall <- 
    map_dfr(all_traits_list, 
          find_evidence_ratios, 
          "R_overall > 0") %>% 
  rename(Overall_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_diff <- 
    map_dfr(all_traits_list, 
          find_evidence_ratios, 
          "R_diff > 0") %>% 
  rename(Diff_Evid_Ratio = Evid.Ratio)
}

Next, we manually calculate evidence ratios for sexually concordant responses to selection by:

1. Using the p_direction function from the bayestestR package, we find the proportion of the posterior distribution (the posterior probability) that is of the median’s sign for each trait in each sex.

2. Using the p_direction output, we calculate the probability that a trait has positive \(R\), \(P(pos)\) or negative \(R\), \(P(neg) = 1 - P(pos)\).

3. Find \(P(concord) = P(pos)_f P(pos)_m + P(neg)_f P(neg)_m\)

4. Find \(P(antag) = P(pos)_f P(neg)_m + P(neg)_f P(pos)_m\)

5. Take the ratio of \(P(concord)\) and \(P(antag)\)

Evidence ratios > 1 indicate that the response to selection is more likely to be sexually concordant, whereas ratios < 1 indicate a sexually antagonistic response has higher probability. As indicated by step 3, a concordant response is possible when the trait responds to selection in either a negative or positive direction in both sexes. Similarly, step 4 shows that an antagonistic response is also possible via two paths. Hence, we once again consider evidence ratios > 39 or < 1/39 as biologically notable (without accounting for multiple testing).

pd_function <-function(Trait_name){
  R_all_traits %>% filter(Trait == Trait_name) %>%
    select(Trait, Response_to_selection_female, Response_to_selection_male) %>%
    p_direction(Response_to_selection_female, method = "direct", null = 0) %>%
    as_tibble() %>%
    pivot_wider(names_from = Parameter, values_from = pd) %>%
    mutate(Trait = Trait_name) %>%
    rename(pd_female = Response_to_selection_female,
           pd_male = Response_to_selection_male)
}

if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){

pd_data <- map_dfr(all_traits_list, pd_function)

Trait_medians <-
  R_all_traits %>%
  group_by(Trait) %>%
  summarise(median_female = median(Response_to_selection_female),
            median_male = median(Response_to_selection_male)) %>%
  ungroup()

evidence_ratios_concord <-
  left_join(pd_data, Trait_medians) %>%
  mutate(prob_pos_female = if_else(median_female > 0, pd_female, 1 - pd_female),
         prob_pos_male = if_else(median_male > 0, pd_male, 1 - pd_male)) %>%
  select(Trait, prob_pos_female, prob_pos_male) %>%

    # Calculate the probabilities that beta_i and beta_j have the same/opposite signs
    mutate(p_sex_concord = prob_pos_female  * prob_pos_male +
             (1 - prob_pos_female) * (1 - prob_pos_male),
           p_sex_antag = prob_pos_female * (1 - prob_pos_male) +
             (1 - prob_pos_female) * prob_pos_male) %>%
    # Find the ratios of these two probabilities (i.e. the "evidence ratios")
    mutate(Concord_Evid_Ratio = p_sex_concord / p_sex_antag)
}

Join the evidence ratios into a single tibble and save as a .csv file for fast loading. Create new columns that indicate whether the evidence ratio indicate biologically notability.

if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){
organismal_level_evidence_ratios <-
  R_evidence_ratios_female %>% 
  left_join(R_evidence_ratios_male) %>% 
  left_join(R_evidence_ratios_overall) %>% 
  left_join(R_evidence_ratios_diff) %>% 
  left_join(evidence_ratios_concord %>% select(Trait, Concord_Evid_Ratio)) %>% 
  mutate(Female_notable = if_else(Female_Evid_Ratio > 39 | Female_Evid_Ratio < 1/39, "Yes", "No"),
         Male_notable = if_else(Male_Evid_Ratio > 39 | Male_Evid_Ratio < 1/39, "Yes", "No"),
         Overall_notable = if_else(Overall_Evid_Ratio > 39 | Overall_Evid_Ratio < 1/39, "Yes", "No"),
         Diff_notable = if_else(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39, "Yes", "No"),
         Concord_notable = if_else(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39, "Yes", "No")) %>% 
  left_join(clean_meta_data) %>% 
  left_join(R_all_traits %>% distinct(Trait, Trait_Sex)) %>%
  mutate(Trait = fct_reorder(Trait, `Trait guild`)) %>%
  arrange(`Trait guild`, Trait) %>% 
  mutate(position = 1:n(),
         Trait_Sex = if_else(Trait_Sex == "Female", "Traits measured in females",
                                         "Traits measured in males")) %>% 
  select(Trait, everything())

write_csv(evidence_ratios, "data/organismal_level_output/organismal_level_evidence_ratios.csv")
} else organismal_level_evidence_ratios <- read_csv("data/organismal_level_output/organismal_level_evidence_ratios.csv")

\(~\)

Build the evidence ratio plots

# work out where to place the x axis labels - we want them in the middle of their trait guild

x_labels <- organismal_level_evidence_ratios %>% 
  group_by(`Trait guild`) %>% 
  summarise(position = min(position) + (max(position) - min(position))/2)

# get colours for each trait guild

guild_pal <- c("#6e7cb9", "#7bbcd5", "#d0e2af", "#f5db99",
               "#e89c81", "#d2848d", "#6e7cb9", "#7bbcd5", 
               "#d0e2af","#f5db99","#e89c81", "#d2848d", 
               "#6e7cb9", "#7bbcd5")

plot_a <- 
  organismal_level_evidence_ratios %>% 
  ggplot(aes(x = position, y = log2(Overall_Evid_Ratio))) +
  geom_point(aes(colour = `Trait guild`),
             size = 4.5, alpha = 1, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12), 
                  labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
  labs(x = NULL,
       #y = "R~F~ log2(ER)",
       y = "_R_  log~2~(ER)") +
  geom_label_repel(data = organismal_level_evidence_ratios %>% 
                     filter(Overall_Evid_Ratio > 200 | Overall_Evid_Ratio < 1/200),  
                   aes(label = Trait_nice),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))
    
plot_b <- 
  organismal_level_evidence_ratios %>% 
  ggplot(aes(x = position, y = log2(Female_Evid_Ratio))) +
  geom_point(aes(colour = `Trait guild`),
             size = 4.5, alpha = 1, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12), 
                  labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
  labs(x = NULL,
       #y = "R~F~ log2(ER)",
       y = "_R_ on females log~2~(ER)") +
  geom_label_repel(data = organismal_level_evidence_ratios %>% 
                     filter(Female_Evid_Ratio > 100 | Female_Evid_Ratio < 1/100),  
                   aes(label = Trait_nice),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))

plot_c <- 
  organismal_level_evidence_ratios %>% 
  ggplot(aes(x = position, y = log2(Male_Evid_Ratio))) +
  geom_point(aes(colour = `Trait guild`),
             size = 4.5, alpha = 1, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12), 
                  labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
  labs(x = "Traits",
       y = "_R_ on males log~2~(ER)") +
  geom_label_repel(data = organismal_level_evidence_ratios %>% 
                     filter(Male_Evid_Ratio > 100 | Male_Evid_Ratio < 1/100),
                   aes(label = Trait_nice),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_text(angle = 90, 
                                 vjust = 0.5,
                                 size = 10),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))

(Figure_S3 <- plot_a / plot_b / plot_c +
  plot_annotation(tag_levels = 'a'))

Past versions of unnamed-chunk-13-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11
b682a4f	ausevo	2023-03-03
6d630f0	ausevo	2022-05-12
7bf0ee0	ausevo	2022-05-10
043902f	tkeaney	2021-07-28
83a0c4f	tkeaney	2021-07-28
32f107f	tkeaney	2021-07-23

Figure S3. Triangles show evidence ratios (ERs) for a overall \(R\), b \(R\) on females (\(R_{FF}\) and \(R_{MF}\)) and c \(R\) on males (\(R_{FM}\) and \(R_{MM}\)). Traits are loosely organised into guilds to aid visualisation. The dashed lines indicate an evidence ratio of 1, where the probability that \(R\) > 0 is equal to the probability that \(R\) < 0. ER > 1 indicates a positive response is more likely, whereas ER < 1 indicates a negative response is more likely. The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.

plot_d <- 
  organismal_level_evidence_ratios %>% 
  ggplot(aes(x = position, y = log2(Diff_Evid_Ratio))) +
  geom_point(aes(colour = `Trait guild`),
             size = 4.5, alpha = 1, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10), 
                     labels = c(paste("1/",2 ^ c(10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10)), limits = c(-12, 12)) +
  labs(x = NULL,
       y = "&Delta;_R_  log~2~(ER)") +
  geom_label_repel(data = organismal_level_evidence_ratios %>% 
                     mutate(Trait_nice = case_when(Trait_nice == "DDT resistance mortality" ~ "DDT susceptibility", .default = Trait_nice)) %>% 
                     filter(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39),  
                   aes(label = Trait_nice),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))

plot_e <- 
  organismal_level_evidence_ratios %>% 
  ggplot(aes(x = position, y = log2(Concord_Evid_Ratio))) +
  geom_point(aes(colour = `Trait guild`),
             size = 4.5, alpha = 1, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10), 
                  labels = c(paste("1/",2 ^ c(10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10)), limits = c(-12, 12)) +
  labs(x = "Trait",
       y = "Sexual concordance log~2~(ER)") +
  geom_label_repel(data = organismal_level_evidence_ratios %>% 
                     filter(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39),  
                   aes(label = Trait_nice),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_text(angle = 90, 
                                 vjust = 0.5,
                                 size = 10),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))

(Figure_3 <- plot_d/ plot_e + plot_annotation(tag_levels = "a"))

Past versions of unnamed-chunk-14-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11
043902f	tkeaney	2021-07-28
83a0c4f	tkeaney	2021-07-28

Figure 3. Triangles show evidence ratios (ERs) for a \(\Delta R\) (\(R\) on females - \(R\) on males) and b sexual concordance, for all measured traits in the organismal level phenotype dataset. Traits are loosely organised into guilds to aid visualisation. The dashed lines indicate an evidence ratio of 1, where the probability that \(\Delta R\) > 0 is equal to the probability that \(\Delta R\) < 0, or that a trait is just as likely to have a sexually concordant response to selection (ER > 1) as a sexually antagonistic response (ER < 1). The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.

Table S2

Table S2 / Supplementary dataset S1: \(R\) partitions and associated evidence ratios for the organismal level trait dataset. Trait descriptions are provided in Table S1. Data can be downloaded as a .csv file from the github repository associated with this project.

all_ol_traits_summary <-
  R_all_traits %>% 
  group_by(Trait, Trait_Sex) %>% 
  median_qi(Response_to_selection_female,
            Response_to_selection_male,
            R_overall,
            R_diff) %>% 
  rename(`Sex trait was measured in` = Trait_Sex,
         R = R_overall,
         `R 2.5% CI` = R_overall.lower,
         `R 97.5% CI` = R_overall.upper,
         `R female` = Response_to_selection_female,
         `R female 2.5% CI` = Response_to_selection_female.lower,
         `R female 97.5% CI` = Response_to_selection_female.upper,
         `R male` = Response_to_selection_male,
         `R male 2.5% CI` = Response_to_selection_male.lower,
         `R male 97.5% CI` = Response_to_selection_male.upper,
         `Delta R` = R_diff,
         `Delta R 2.5% CI` = R_diff.lower,
         `Delta R 97.5% CI` = R_diff.upper) %>% 
  left_join(organismal_level_evidence_ratios %>% select(-contains("notable")) %>% 
              rename(`R ER` = Overall_Evid_Ratio, `R female ER` = Female_Evid_Ratio,
                     `R male ER` = Male_Evid_Ratio, `Delta R ER` = Diff_Evid_Ratio,
                     `Sexual concordance ER` = Concord_Evid_Ratio)) %>% 
  select(Trait, Trait_nice, `Sex trait was measured in`, R, `R 2.5% CI`, `R 97.5% CI`, `R ER`, 
         `R female`, `R female 2.5% CI`, `R female 97.5% CI`, `R female ER`,
         `R male`, `R male 2.5% CI`, `R male 97.5% CI`, `R male ER`,
         `Delta R`, `Delta R 2.5% CI`, `Delta R 97.5% CI`, `Delta R ER`, `Sexual concordance ER`) %>%
  ungroup() %>% 
  mutate(across(c(4,5,6,8,9,10,12,13,14,16,17,18), ~round(.x, 3)),
         across(ends_with("ER"), ~ round(.x, 4)))

write_csv(all_ol_traits_summary, "data/organismal_level_output/Supplementary_dataset_1.csv")

my_data_table(all_ol_traits_summary)

\(~\)

Estimate \(\overline{|R|}\) in each sex

The model

Fit the model to test whether \(|R|\) is affected by the sex fitness and trait values were measured in. Here, we take a meta-analysis type approach and weight each \(R\) estimate by the inverse of its standard error.

R_medians <-
  R_long_form %>%
  group_by(Trait, Fitness_Sex, Trait_Sex) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>% 
  ungroup() %>% 
  select(-variable) %>% 
  left_join(clean_meta_data) %>% 
  mutate(absolute_R = abs(median),
         Fitness_Sex = as.factor(Fitness_Sex),
         Trait_Sex = as.factor(Trait_Sex),
         Trait = as.factor(Trait)) 

median_R_model <- 
  brm(absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1|study_ID) + (1|Trait),
      family = brmsfamily(family = "Gamma"), # gamma is appropriate for the half-normal distribution created by taking the absolute
      data = R_medians,
      prior = c(prior(normal(-2.2, 1), class = Intercept),
                prior(exponential(1), class = sd),
                prior(exponential(1), class = shape),
                prior(normal(0, 0.25), class = b)),
      warmup = 4000, iter = 12000,
      seed = 1, cores = 4, chains = 4,
      control = list(adapt_delta = 0.9, max_treedepth = 10),
      file = "fits/median_R_model")

print(median_R_model)

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1 | study_ID) + (1 | Trait) 
   Data: R_medians (Number of observations: 950) 
  Draws: 4 chains, each with iter = 12000; warmup = 4000; thin = 1;
         total post-warmup draws = 32000

Group-Level Effects: 
~study_ID (Number of levels: 76) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.32      0.05     0.23     0.43 1.00     6142    10131

~Trait (Number of levels: 475) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.59      0.02     0.54     0.63 1.00     8421    12837

Population-Level Effects: 
                              Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                        -2.70      0.06    -2.82    -2.58 1.00
Fitness_SexMale                   0.18      0.02     0.14     0.22 1.00
Trait_SexMale                    -0.15      0.06    -0.27    -0.03 1.00
Fitness_SexMale:Trait_SexMale     0.11      0.03     0.05     0.17 1.00
                              Bulk_ESS Tail_ESS
Intercept                        16404    18313
Fitness_SexMale                  36862    23998
Trait_SexMale                    11981    15341
Fitness_SexMale:Trait_SexMale    34545    23661

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     2.42      0.03     2.36     2.49 1.00    39872    21791

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

To check how does the gamma distribution fares at capturing the shape of the data, we use the model to recapitulate the data it was trained on. Given the shape of the predicted data is similar to the real data, there’s a high chance we’ve used an appropriate distribution family.

pp_check(median_R_model)

Past versions of unnamed-chunk-17-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11

\(~\)

Build Panels a and b for Figure 2

Get model predictions and plot

new_data <- expand_grid(Fitness_Sex = c("Female", "Male"),
                        Trait_Sex = c("Female", "Male"))

R_fitted <-
  fitted(median_R_model, newdata = new_data, re_formula = NA, summary = F) %>% 
  as.data.frame() %>% 
  rename(FemaleFitness_FemaleTrait = V1, FemaleFitness_MaleTrait = V2, 
         MaleFitness_FemaleTrait = V3, MaleFitness_MaleTrait = V4) %>% 
  as_tibble() %>% 
  mutate(percent_diff_female = ((MaleFitness_FemaleTrait - FemaleFitness_FemaleTrait) / FemaleFitness_FemaleTrait)*100,
         percent_diff_male = ((MaleFitness_MaleTrait - FemaleFitness_MaleTrait)/ FemaleFitness_MaleTrait)*100) %>% 
  pivot_longer(cols = everything(), names_to = "Parameter", values_to = "R_mean") %>% 
  mutate(Fitness_Sex = case_when(str_detect(Parameter, "FemaleFitness") ~ "Female",
                                 str_detect(Parameter, "MaleFitness") ~ "Male"),
         Trait_Sex = case_when(str_detect(Parameter, "FemaleTrait") ~ "Female",
                               str_detect(Parameter, "MaleTrait") ~ "Male",
                               str_detect(Parameter, "percent_diff_female") ~ "Female",
                               str_detect(Parameter, "percent_diff_male") ~ "Male"))

R_p1 <-
  R_fitted %>% 
  filter(!str_detect(Parameter, "percent")) %>% 
  ggplot(aes(x = R_mean, y = Trait_Sex, fill = Fitness_Sex)) + #fill = Fitness_Sex)) +
  stat_slab(alpha = 0.9, shape = 21) +#,
  labs(x = expression("Mean |"* italic(R) * "| for organismal level traits"),
       y = "Phenotyped sex") +
  scale_fill_manual(values = c(carto_pal(7, "Purp")[5], carto_pal(7, "Peach")[5])) +
  coord_cartesian(xlim = c(0.045, 0.10)) +
  theme_minimal() + 
  theme(panel.background = element_rect(fill='transparent'),
        #panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA),
        legend.position = "none",
        text = element_text(size=13))

R_p2 <-
  R_fitted %>% 
  filter(str_detect(Parameter, "percent")) %>%
  ggplot(aes(x = R_mean, y = Trait_Sex, fill = after_stat(x > 0))) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, fill = carto_pal(7, "Emrld")[2],
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.2) + 
  geom_vline(xintercept = 0, linetype = 2, linewidth = 1.2) +
  #coord_cartesian(xlim = c(-5, 50)) +
  xlab(expression("% sex difference in |"*italic(R)* "|")) +
  ylab("Phenotyped sex") +
  coord_cartesian(xlim = c(-5, 45)) +
  theme_minimal() + 
  theme(panel.background = element_rect(fill='transparent'),
        #panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA),
        legend.position = "none",
        text = element_text(size=13))

R_p1 + R_p2 +
  plot_annotation(tag_levels = 'a')

Past versions of unnamed-chunk-18-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11
daa935d	ausevo	2023-03-08
b682a4f	ausevo	2023-03-03

\(~\)

Calculating stats for the results section

R_fitted %>% 
  group_by(Parameter) %>% 
  median_qi(R_mean)

# A tibble: 6 × 7
  Parameter                  R_mean  .lower  .upper .width .point .interval
  <chr>                       <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 FemaleFitness_FemaleTrait  0.0673  0.0598  0.0761   0.95 median qi       
2 FemaleFitness_MaleTrait    0.0579  0.0513  0.0653   0.95 median qi       
3 MaleFitness_FemaleTrait    0.0808  0.0716  0.0915   0.95 median qi       
4 MaleFitness_MaleTrait      0.0776  0.0689  0.0876   0.95 median qi       
5 percent_diff_female       20.1    15.3    25.1      0.95 median qi       
6 percent_diff_male         34.1    28.7    39.7      0.95 median qi       

\(~\)

Transcriptome analysis

Selecting data appropriate for analysis

We use the data archived on the dgrp2 website, which was provided by Huang et al. (2015). These data are levels of expression of genes measured across 185 DGRP lines. Huang et al’s data contains Y-linked genes that have higher/equal expression in females in all lines, presumably microarray issues/errors. To be conservative, we restrict our analyses to genes that are known to be on chromosomes that are present in both sexes. After data cleaning, we retain 14,268 genes in our analysis.

We also load gene annotation data using the org.Dm.eg.db R package provided by BiocManager. The code used to produce annotations is provided in the code subdirectory, in the get_gene_annotations.R file.

Huang et al. 2015 replicated each gene expression measurement twice. To find line means, we simply take the average value for each line. We then standardise the expression of each gene to have \(\mu = 0\) and \(\sigma = 1\) (as with the ‘phenotypic’ traits.

# load in the gene names, we'll use these for plots and tables

Dmel_gene_names <- read_csv("data/input/all_dmel_genes.csv")


gene_annotations <- read_csv("data/input/gene_anntotations.csv")


# Helper to load all the Huang et al. expression data into tidy format
load_expression_data <- function(gene_annotations){
  
  # Note: Huang et al's data contains Y-linked genes that have
  # higher/equal expression in *females* in all lines, presumably microarray issues/errors.
  # To be conservative, we restrict our analyses to genes that are known to be on a
  # chromosomes that is present in both sexes
  genes_allowed <- gene_annotations %>%
    filter(chromosome %in% c("2L", "2R", "3L", "3R", "4", "X")) %>%
    pull(FBID) 
  
  females <- read_delim("data/input/huang_transcriptome/dgrp.array.exp.female.txt", delim = " ") %>%
    filter(gene %in% genes_allowed) %>% 
    gather(sample, expression, -gene) %>% 
    mutate(line = map_chr(str_extract_all(sample, "line_[:digit:]*"), ~ .x[1]),
           replicate = map_chr(str_split(sample, ":"), ~ .x[2]),
           sex = "Female")
  
  males <- read_delim("data/input/huang_transcriptome/dgrp.array.exp.male.txt", delim = " ") %>%
    filter(gene %in% genes_allowed) %>% 
    gather(sample, expression, -gene) %>% 
    mutate(line = map_chr(str_extract_all(sample, "line_[:digit:]*"), ~ .x[1]),
           replicate = map_chr(str_split(sample, ":"), ~ .x[2]),
           sex = "Male")
  
  bind_rows(females, males) %>% 
    select(sample, line, sex, replicate, gene, expression)
}

expression_line_means <- 
  load_expression_data(gene_annotations) %>% # use the custom function to load expression data
  mutate(line = str_remove(line, "line_"),
         line = as.factor(line)) %>% 
  group_by(line, sex, gene) %>% 
  mutate(trait_value = mean(expression)) %>% # compute the average between replicates for each gene
  ungroup() %>% 
  distinct(line, sex, gene, trait_value) %>% 
  group_by(sex, gene) %>% # scale the traits, specific to gene and sex
  mutate(trait_value = as.numeric(scale(trait_value))) %>% 
  ungroup() %>% 
  rename(Sex = sex) %>% 
# join the fitness data
left_join(trait_data %>% distinct(line, female_fitness, male_fitness))

# the transcriptome is large; memory is therefore a consistent problem with this analysis - it helps to clear often
# 
invisible(gc())

\(~\)

Calculating R: the response to selection

\(~\)

Build models to calculate \(R_{FF}\) & \(R_{MF}\)

To estimate the covariance between \(A_w\) and \(A_z\) (which in this case is a transcript), we fitted bivariate Bayesian linear models using the brms package (Bürkner, 2017) for R version 4.2.2. For each combination of trait/transcript and sex, we used line means for the focal trait/transcript and the fitness of the focal sex as the two response variables and fitted an intercept-only Gaussian model. Each model returned a posterior distribution of the residual correlation between trait/transcript and fitness, which for data expressed in standard units is equivalent to the covariance.

Build functions to run the models

# RFF estimates

female_transcripts <- expression_line_means %>% filter(Sex == "Female")

transcript_list <- unique(female_transcripts$gene) # an input to the map_dfr() function that we'll need in a few chunks time

# code the model structure we will use for all traits/transcripts using one example - `FBgn0000014`. We can then use the update() function to run this model many times, once for each trait/transcript measured in females. update() makes this process many times faster, because the model can immediately start sampling, without the need to recompile. 

RFF_transcriptome_model <-
  brm(data = female_transcripts %>% filter(gene == "FBgn0000014"),
      family = gaussian,
      bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = femalefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, file = "fits/RFF_transcriptome_model",
      backend = "cmdstanr", refresh = 400)  


# make a function to update the model and the posterior sample output with the 'selected trait'

RFF_transcriptome_calculator <- function(selected_gene){
  
  data <- female_transcripts %>% filter(gene == selected_gene)
  
  model <- update(
    RFF_transcriptome_model, newdata = data,
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, backend = "cmdstanr", refresh = 400)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>% 
    mutate(Trait = selected_gene) %>% 
    select(Trait, Response_to_selection_female) %>% 
    as_tibble()
  
  posterior
}

# RMF estimates

RMF_transcriptome_model <-
  brm(data = female_transcripts %>% filter(gene == "FBgn0000014"),
      family = gaussian,
      bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = malefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, file = "fits/RMF_transcriptome_model",
      backend = "cmdstanr", refresh = 400) 

# make a function to update the model and the posterior sample output with your desired trait

RMF_transcriptome_calculator <- function(selected_gene){
  
  data <- female_transcripts %>% filter(gene == selected_gene)
  
  model <- update(
    RMF_transcriptome_model, newdata = data,
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, backend = "cmdstanr", refresh = 400)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>% 
    mutate(Trait = selected_gene) %>% 
    select(Trait, Response_to_selection_male) %>% 
    as_tibble()
  
  posterior
}

\(~\)

Build models to calculate \(R_{FM}\) & \(R_{MF}\)

\(~\)

# RMM estimates

male_transcripts <- expression_line_means %>% filter(Sex == "Male")

RMM_transcriptome_model <-
  brm(data = male_transcripts %>% filter(gene == "FBgn0000014"),
      family = gaussian,
      bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = malefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, file = "fits/RMM_transcriptome_model",
      backend = "cmdstanr", refresh = 400)    

# make a function to update the model and the posterior sample output with your desired trait

RMM_transcriptome_calculator <- function(selected_gene){
  
  data <- male_transcripts %>% filter(gene == selected_gene)
  
  model <- update(
    RMM_transcriptome_model, newdata = data,
    chains = 4, cores = 4, iter = 4000, warmup = 2000,
    seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
    mutate(Trait = selected_gene) %>% 
    select(Trait, Response_to_selection_male) %>% 
    as_tibble()
  
  posterior
}

# RFM estimates

RFM_transcriptome_model <-
  brm(data = male_transcripts %>% filter(gene == "FBgn0000014"),
      family = gaussian,
      bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
      prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
                prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
                prior(normal(1, 0.1), class = sigma, resp = femalefitness),
                prior(normal(1, 0.1), class = sigma, resp = traitvalue),
                prior(lkj(2), class = rescor)),
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1, file = "fits/RFM_transcriptome_model",
      backend = "cmdstanr", refresh = 400)    

# make a function to update the model and the posterior sample output with your desired trait

RFM_transcriptome_calculator <- function(selected_trait){
  
  data <- male_transcripts %>% filter(gene == selected_trait)
  
  model <- update(
    RFM_transcriptome_model, newdata = data,
      chains = 4, cores = 4, iter = 4000, warmup = 2000,
      seed = 1)
  
  posterior <- 
    as_draws_df(model) %>% 
    rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>% 
    mutate(Trait = selected_gene) %>% 
    select(Trait, Response_to_selection_female) %>% 
    as_tibble()
  
  posterior
}

\(~\)

Run the models for all the traits

Run the models using RFF_transcriptome_calculator, RMF_transcriptome_calculator, RMM_transcriptome_calculator and RFM_transcriptome_calculator.

This takes a fair bit of memory, so it might be necessary to run the models in chunks rather than all in one go. To do this, you can break the transcript_list_female list into parts and feed it into the map_dfr function. The completed subset can then be saved to your hard disk, removed from R and cleared from your computers memory. This frees up memory to run another chunk without losing progress. Expand the code chunk below to see an example. All other chunks have been run and saved for later use.

# run the RFF function

transcript_list_1 <- transcript_list[1:2000]

if(!file.exists("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")){
  RFF_transcript_1 <- map_dfr(transcript_list_1, RFF_transcriptome_calculator)
  write_csv(RFF_transcript_1, 
            file = "data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")
  rm(RFF_transcript_1)
  invisible(gc())
} else RFF_transcript_1 <- read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")

Load all the posterior draws, combine and summarise and save the result to the hard disk. This allows us to simply load the summarised results once everything has been run once.

if(!file.exists("data/transcriptome_output/R_summarised_transcriptome.csv")){
RFF_transcriptome_complete <- 
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv")
  ) %>% 
  group_by(Trait) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
  ungroup() %>% 
  select(-variable) %>% 
  mutate(absolute_R = abs(median),
         Fitness_Sex = "Female",
         Trait_Sex = "Female")

invisible(gc())

RMF_transcriptome_complete <- 
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv")
  ) %>% 
  group_by(Trait) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
  ungroup() %>% 
  select(-variable) %>% 
  mutate(absolute_R = abs(median),
         Fitness_Sex = "Male",
         Trait_Sex = "Female")

invisible(gc())

RMM_transcriptome_complete <- 
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv")
  ) %>% 
  group_by(Trait) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
  ungroup() %>% 
  select(-variable) %>% 
  mutate(absolute_R = abs(median),
         Fitness_Sex = "Male",
         Trait_Sex = "Male")

invisible(gc())

RFM_transcriptome_complete <- 
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv")
  ) %>% 
  group_by(Trait) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
  ungroup() %>% 
  select(-variable) %>% 
  mutate(absolute_R = abs(median),
         Fitness_Sex = "Female",
         Trait_Sex = "Male")

invisible(gc())

R_summarised_transcriptome <-
  bind_rows(RFF_transcriptome_complete, RMF_transcriptome_complete,
            RFM_transcriptome_complete, RMM_transcriptome_complete)

write_csv(R_summarised_transcriptome, "data/transcriptome_output/R_summarised_transcriptome.csv")

} else R_summarised_transcriptome <- read_csv("data/transcriptome_output/R_summarised_transcriptome.csv")

\(~\)

Calculate various metrics

First, find \(R_{overall}\)

if(!file.exists("data/transcriptome_output/R_overall_transcript_summarised.csv")){
  R_overall_transcript_summarised <-
    bind_rows(
      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv(
                  "data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      # now the traits measured in males
      
      left_join(read_csv(
        "data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
        select(Trait, R_overall) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_overall = median) %>%
        select(-variable) %>% 
        mutate(Trait_Sex = "Male")
    )
  
  write_csv(R_overall_transcript_summarised, "data/transcriptome_output/R_overall_transcript_summarised.csv")
} else R_overall_transcript_summarised <- read_csv("data/transcriptome_output/R_overall_transcript_summarised.csv")

Now find \(\Delta R\)

if(!file.exists("data/transcriptome_output/R_diff_transcript_summarised.csv")){
  R_diff_transcript_summarised <-
    bind_rows(
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
                mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
 
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>%  
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),

      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Female"),
      
      # now the traits measured in males
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male"),
      
      left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(),
                read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>% 
                  group_by(Trait) %>% 
                  mutate(draw = row_number()) %>% 
                  ungroup(), by = c("Trait", "draw")) %>% 
        mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
        select(Trait, R_diff) %>% 
        group_by(Trait) %>% 
        summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
        ungroup() %>% 
        rename(R_diff = median) %>% 
        select(-variable) %>% 
        mutate(Trait_Sex = "Male")
    )
  
  write_csv(R_diff_transcript_summarised, "data/transcriptome_output/R_diff_transcript_summarised.csv")
} else R_diff_transcript_summarised <- read_csv("data/transcriptome_output/R_diff_transcript_summarised.csv")

Calculating evidence ratios

Use the brms hypothesis() function to compute evidence ratios, just as we did for the organismal level phenotypic traits.

# build a function to calculate the evidence ratio, get the relevant info from the output and convert it to a tibble. Basically the same function, just with a different dataframe input and a memory clearing step that slows it down but allows it to run for many traits. 

find_evidence_ratios_transcriptome <- function(Trait_name, specify_hypothesis){
  x <- hypothesis(R_transcriptome %>% 
                    filter(Trait == Trait_name),
                  specify_hypothesis)
  
  x <- x$hypothesis
  
  invisible(gc())
  
  x %>% as_tibble() %>% 
    select(Evid.Ratio) %>% 
    mutate(Trait = Trait_name)
}

This pc has 16gb memory; not enough to load in all the transcript posterior draws at once. Having previously saved these to the hard disk we load the draws for all female measured traits and calculate the evidence ratio that \(R_{FF} > 0\), \(R_{MF} > 0\), \(R_F > 0\) and \(\Delta R > 0\). The draws for the female measured traits were then removed from the computers memory using the rm() and gc() functions and the process was repeated for the male measured traits.

# find evidence ratios for the traits measured in females

if(!file.exists("data/transcriptome_output/evidence_ratios_transcriptome_female.csv")){ # name refers to traits measured in females

R_transcriptome <-
  cbind(
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv")
  ),
   rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv")) %>% 
      select(-Trait)
   )

R_transcriptome <- 
  R_transcriptome %>% 
  mutate(R_overall = (Response_to_selection_female - Response_to_selection_male) / 2,
         R_diff = Response_to_selection_female - Response_to_selection_male) 
  
# calculate evidence ratios 

R_evidence_ratios_female <- 
  map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "Response_to_selection_female > 0") %>% 
  rename(Female_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_male <-
  map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "Response_to_selection_male > 0") %>% 
  rename(Male_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_overall <- 
    map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "R_overall > 0") %>% 
  rename(Overall_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_diff <- 
    map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "R_diff > 0") %>% 
  rename(Diff_Evid_Ratio = Evid.Ratio)

# write each result to csv in 

evidence_ratios_transcriptome_female <-
  R_evidence_ratios_female %>% 
  left_join(R_evidence_ratios_male) %>% 
  left_join(R_evidence_ratios_overall) %>% 
  left_join(R_evidence_ratios_diff) %>% 
  select(Trait, everything()) %>% 
  rename(Female_Evid_Ratio = Female_evidence_ratio) %>%  # fixing small typo
  mutate(Trait_Sex = "Female")

write_csv(evidence_ratios_transcriptome_female, "data/transcriptome_output/evidence_ratios_transcriptome_female.csv")

} else evidence_ratios_transcriptome_female <- read_csv("data/transcriptome_output/evidence_ratios_transcriptome_female.csv")

# find evidence ratios for the traits measured in males

if(!file.exists("data/transcriptome_output/evidence_ratios_transcriptome_male.csv")){ # name refers to traits measured in females

R_transcriptome <-
  cbind(
  rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv")
  ),
   rbind(
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv"),
    read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv")) %>% 
      select(-Trait)
   )

R_transcriptome <- 
  R_transcriptome %>% 
  mutate(R_overall = (Response_to_selection_female - Response_to_selection_male) / 2,
         R_diff = Response_to_selection_female - Response_to_selection_male) 
  
# calculate evidence ratios 

R_evidence_ratios_female <- 
  map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "Response_to_selection_female > 0") %>% 
  rename(Female_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_male <-
  map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "Response_to_selection_male > 0") %>% 
  rename(Male_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_overall <- 
    map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "R_overall > 0") %>% 
  rename(Overall_Evid_Ratio = Evid.Ratio)

R_evidence_ratios_diff <- 
    map_dfr(transcript_list, 
          find_evidence_ratios_transcriptome, 
          "R_diff > 0") %>% 
  rename(Diff_Evid_Ratio = Evid.Ratio)

evidence_ratios_transcriptome_male <-
  R_evidence_ratios_female %>% 
  left_join(R_evidence_ratios_male) %>% 
  left_join(R_evidence_ratios_overall) %>% 
  left_join(R_evidence_ratios_diff) %>% 
  select(Trait, everything()) %>% 
  rename(Female_Evid_Ratio = Female_evidence_ratio) %>%  # fixing small typo
  mutate(Trait_Sex = "Male")

write_csv(evidence_ratios_transcriptome_male, "data/transcriptome_output/evidence_ratios_transcriptome_male.csv")

} else evidence_ratios_transcriptome_male <- read_csv("data/transcriptome_output/evidence_ratios_transcriptome_male.csv")

To calculate evidence ratios for sexually concordant responses to selection we use the same function as for the organismal-level phenotypic traits. Run the function in chunks, save them to hard-disk, clear memory and repeat. If already done, just load.

if(!file.exists("data/transcriptome_output/concord_evidence_ratio_data.csv")){

female_transcripts_1 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_1 <- unique(female_transcripts_1$Trait)

f1 <- get_evidence_ratios(female_transcripts_1, trait_list_1, "Female")

rm(female_transcripts_1)
invisible(gc())

female_transcripts_2 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_2 <- unique(female_transcripts_2$Trait)

f2 <- get_evidence_ratios(female_transcripts_2, trait_list_2, "Female")

rm(female_transcripts_2)
invisible(gc())

female_transcripts_3 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_3 <- unique(female_transcripts_3$Trait)

f3 <- get_evidence_ratios(female_transcripts_3, trait_list_3, "Female")

rm(female_transcripts_3)
invisible(gc())

female_transcripts_4 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_4 <- unique(female_transcripts_4$Trait)

f4 <- get_evidence_ratios(female_transcripts_4, trait_list_4, "Female")

rm(female_transcripts_4)
invisible(gc())

female_transcripts_5 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_5 <- unique(female_transcripts_5$Trait)

f5 <- get_evidence_ratios(female_transcripts_5, trait_list_5, "Female")

rm(female_transcripts_5)
invisible(gc())

female_transcripts_6 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_6 <- unique(female_transcripts_6$Trait)

f6 <- get_evidence_ratios(female_transcripts_6, trait_list_6, "Female")

rm(female_transcripts_6)
invisible(gc())

female_transcripts_7 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

trait_list_7 <- unique(female_transcripts_7$Trait)

f7 <- get_evidence_ratios(female_transcripts_7, trait_list_7, "Female")

rm(female_transcripts_7)
invisible(gc())

# now the transcripts measured in males

male_transcripts_1 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_1 <- unique(male_transcripts_1$Trait)

m1 <- get_evidence_ratios(male_transcripts_1, m_trait_list_1, "Male")

rm(male_transcripts_1)
invisible(gc())

male_transcripts_2 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_2 <- unique(male_transcripts_2$Trait)

m2 <- get_evidence_ratios(male_transcripts_2, m_trait_list_2, "Male")

rm(male_transcripts_2)
invisible(gc())

male_transcripts_3 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_3 <- unique(male_transcripts_3$Trait)

m3 <- get_evidence_ratios(male_transcripts_3, m_trait_list_3, "Male")

rm(male_transcripts_3)
invisible(gc())

male_transcripts_4 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_4 <- unique(male_transcripts_4$Trait)

m4 <- get_evidence_ratios(male_transcripts_4, m_trait_list_4, "Male")

rm(male_transcripts_4)
invisible(gc())

male_transcripts_5 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_5 <- unique(male_transcripts_5$Trait)

m5 <- get_evidence_ratios(male_transcripts_5, m_trait_list_5, "Male")

rm(male_transcripts_5)
invisible(gc())

male_transcripts_6 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_6 <- unique(male_transcripts_6$Trait)

m6 <- get_evidence_ratios(male_transcripts_6, m_trait_list_6, "Male")

rm(male_transcripts_6)
invisible(gc())

male_transcripts_7 <-
  left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(),
            read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>% 
              group_by(Trait) %>% 
              mutate(draw = row_number()) %>% 
              ungroup(), by = c("Trait", "draw"))

m_trait_list_7 <- unique(male_transcripts_7$Trait)

m7 <- get_evidence_ratios(male_transcripts_7, m_trait_list_7, "Male")

rm(male_transcripts_7)
invisible(gc())

concord_evidence_ratio_data <- 
  bind_rows(f1, f2, f3, f4, f5, f6, f7, m1, m2, m3, m4, m5, m6, m7)

write_csv(concord_evidence_ratio_data, "data/transcriptome_output/concord_evidence_ratio_data.csv")

} else concord_evidence_ratio_data <- read_csv("data/transcriptome_output/concord_evidence_ratio_data.csv")

Bind all the evidence ratios together

evidence_ratios_transcriptome <-
  evidence_ratios_transcriptome_female %>% 
  bind_rows(evidence_ratios_transcriptome_male %>% rename(Male_Evid_Ratio = Male_evidence_ratio)) %>%  # fix a typo
  left_join(concord_evidence_ratio_data %>% rename(Concord_Evid_Ratio = evidence_ratio_concord) %>% select(Trait, Trait_Sex, Concord_Evid_Ratio),
            by = c("Trait", "Trait_Sex")) %>% 
  left_join(gene_annotations %>% rename(Trait = FBID)) %>% 
  mutate(chromosome_numeric = case_when(chromosome == "2L" ~ 1,
                                        chromosome == "2R" ~ 2,
                                        chromosome == "3L" ~ 3,
                                        chromosome == "3R" ~ 4,
                                        chromosome == "X" ~ 5)) %>% 
  #left_join(R_summarised_transcriptome %>% distinct(Trait, Trait_Sex)) %>%
  #mutate(Trait = fct_reorder(Trait, chromosome_numeric))
  arrange(chromosome, Trait) %>% 
  mutate(position = 1:n(),
         Trait_Sex = if_else(Trait_Sex == "Female", "Traits measured in females",
                                         "Traits measured in males")) %>% 
  select(Trait, everything())

\(~\)

Find the number of evidence ratios > 39 or < 1/39 for each partition of \(R\)

response_pos <- evidence_ratios_transcriptome %>% filter(Overall_Evid_Ratio > 39) %>% arrange(-Overall_Evid_Ratio)

response_neg <- evidence_ratios_transcriptome %>% filter(Overall_Evid_Ratio < 1/39) %>% arrange(Overall_Evid_Ratio)

response_to_female_pos <- evidence_ratios_transcriptome %>% filter(Female_Evid_Ratio > 39) %>% arrange(-Female_Evid_Ratio)

response_to_female_neg <- evidence_ratios_transcriptome %>% filter(Female_Evid_Ratio < 1/39) %>% arrange(Female_Evid_Ratio)

response_to_male_pos <- evidence_ratios_transcriptome %>% filter(Male_Evid_Ratio > 39) %>% arrange(-Male_Evid_Ratio)

response_to_male_neg <- evidence_ratios_transcriptome %>% filter(Male_Evid_Ratio < 1/39) %>% arrange(Male_Evid_Ratio)

response_sexes_diff <- 
  evidence_ratios_transcriptome %>% filter(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39) %>% arrange(-Diff_Evid_Ratio)

sexually_concordant <- evidence_ratios_transcriptome %>% filter(Concord_Evid_Ratio > 39) %>% arrange(-Concord_Evid_Ratio)

sexually_antagonistic <- evidence_ratios_transcriptome %>% filter(Concord_Evid_Ratio < 1/39) %>% arrange(Concord_Evid_Ratio)

684 transcripts are predicted to have a positive overall response to selection, while 844 transcripts ar expected to have a negative overall response.

In response to selection on females, 424 transcripts are expected to increase in expression, whereas 530 are expected to decrease in mean expression.

In response to selection on males, 468 transcripts are expected to increase in expression, whereas 619 are expected to decrease in mean expression.

575 transcripts have notably different responses to selection on females compared with males.

1 transcripts are expected to have a sexually antagonistic response to selection, while 25 transcripts are expected to have a sexually concordant response. From these numbers, we can roughly estimate that 3.8% of transcripts that are responding to selection in both sexes are sexually antagonistic.

Build the evidence ratio plots

x_labels <- evidence_ratios_transcriptome %>% 
  group_by(chromosome) %>% 
  summarise(position = min(position) + (max(position) - min(position))/2)


guild_pal <- c("#6e7cb9", "#7bbcd5", "#d0e2af", "#f5db99",
               "#e89c81", "#d2848d")

plot_transcriptome_a <- 
  evidence_ratios_transcriptome %>% 
  ggplot(aes(x = position, y = log2(Overall_Evid_Ratio))) +
  geom_point(aes(colour = chromosome),
             size = 2.5, alpha = 0.8, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14), 
                     labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-18, 18)) +
  labs(x = NULL,
       #y = "R~F~ log2(ER)",
       y = "_R_  log~2~(ER)") +
  geom_label_repel(data = evidence_ratios_transcriptome %>% 
                     filter(Overall_Evid_Ratio > 1000 | Overall_Evid_Ratio < 1/1000 ),  
                   aes(label = gene_symbol),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 10))

plot_transcriptome_b <- 
  evidence_ratios_transcriptome %>% 
  ggplot(aes(x = position, y = log2(Female_Evid_Ratio))) +
  geom_point(aes(colour = chromosome),
             size = 2.5, alpha = 0.8, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14), 
                     labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-16, 16)) +
  labs(x = NULL,
       #y = "R~F~ log2(ER)",
       y = "_R_ on females log~2~(ER)") +
  geom_label_repel(data = evidence_ratios_transcriptome %>% 
                     filter(Female_Evid_Ratio > 1000 | Female_Evid_Ratio < 1/1000 ), 
                   aes(label = gene_symbol),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 10))

plot_transcriptome_c <- 
  evidence_ratios_transcriptome %>% 
  ggplot(aes(x = position, y = log2(Male_Evid_Ratio))) +
   geom_point(aes(colour = chromosome),
             size = 2.5, alpha = 0.8, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14), 
                     labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-15, 15)) +
  labs(x = "Chromosome arm",
       #y = "R~F~ log2(ER)",
       y = "_R_ on males log~2~(ER)") +
  geom_label_repel(data = evidence_ratios_transcriptome %>% 
                     filter(Male_Evid_Ratio > 1000 | Male_Evid_Ratio < 1/1000 ),  
                   aes(label = gene_symbol),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_text(size = 15),
        #axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 10))

(Figure_S4 <- plot_transcriptome_a / plot_transcriptome_b / plot_transcriptome_c +
  plot_annotation(tag_levels = 'a')
)

Past versions of unnamed-chunk-33-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11

Figure S4. Triangles show evidence ratios (ERs) for gene expression traits. a shows overall \(R\), b \(R\) on females (\(R_{FF}\) and \(R_{MF}\)) and c \(R\) on males (\(R_{FM}\) and \(R_{MM}\)). The dashed lines indicate an evidence ratio of 1, where the probability that \(R\) > 0 is equal to the probability that \(R\) < 0. ER > 1 indicates a positive response is more likely, whereas ER < 1 indicates a negative response is more likely. The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.

plot_transcriptome_d <- 
  evidence_ratios_transcriptome %>% 
  ggplot(aes(x = position, y = log2(Diff_Evid_Ratio))) +
   geom_point(aes(colour = chromosome),
             size = 2.5, alpha = 0.8, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12), 
                     labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-13, 13)) +
  labs(x = "Chromosome arm",
       #y = "R~F~ log2(ER)",
       y = "&Delta;_R_  log~2~(ER)") +
  geom_label_repel(data = evidence_ratios_transcriptome %>% 
                     filter(Diff_Evid_Ratio > 250 | Diff_Evid_Ratio < 1/250 ),  
                   aes(label = gene_symbol),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_markdown(size = 15))

plot_transcriptome_e <- 
  evidence_ratios_transcriptome %>% 
  ggplot(aes(x = position, y = log2(Concord_Evid_Ratio))) +
   geom_point(aes(colour = chromosome),
             size = 2.5, alpha = 0.8, shape = 17) +
  geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
  geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
  geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
  #scale_shape_manual(values = c(25, 24)) +
  scale_colour_manual(values = guild_pal) +
  scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
  scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12), 
                     labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-12, 12)) +
  labs(x = "Chromosome arm",
       y = "Sexual concrdoance  log~2~(ER)") +
  geom_label_repel(data = evidence_ratios_transcriptome %>% 
                     filter(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39), 
                   aes(label = gene_symbol),
                   fill = "white",
                   size = 3, alpha = 0.9,
                   min.segment.length = 0, seed = 1,
                   box.padding = 0.5, point.padding = 0.5,
                   max.overlaps = 30) +
  theme_minimal() +
    theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        axis.text.x=element_text(size = 15),
        #axis.ticks.x=element_blank(),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_markdown(size = 15))


(Figure_4 <- plot_transcriptome_d/ plot_transcriptome_e + plot_annotation(tag_levels = "a")
)

Past versions of unnamed-chunk-34-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11

Figure 4. Triangles show evidence ratios (ERs) for a \(\Delta R\) (\(R\) on females - \(R\) on males) and b sexual concordance, for all measured traits in the transcriptome dataset. The dashed lines indicate an evidence ratio of 1, where the probability that \(\Delta R\) > 0 is equal to the probability that \(\Delta R\) < 0, or that a trait is just as likely to have a sexually concordant response to selection (ER > 1) as a sexually antagonistic response (ER < 1). The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.

\(~\)

Table S3

Table S3 / Supplementary dataset S2: \(R\) partitions and associated evidence ratios for the gene expression trait dataset. Data can be downloaded as a .csv file from the github repository associated with this project.

all_transcriptome_summary <-
  R_summarised_transcriptome %>% 
  select(-c(sd, absolute_R)) %>% 
  pivot_wider(names_from = c(Fitness_Sex), values_from = 2:4) %>% 
  rename(`R female` = median_Female,
         `R female 2.5% CI` = `2.5%_Female`,
         `R female 97.5% CI` = `97.5%_Female`,
         `R male` = median_Male,
         `R male 2.5% CI` = `2.5%_Male`,
         `R male 97.5% CI` = `97.5%_Male`) %>% 
  left_join(
    
    R_overall_transcript_summarised %>%
      select(-sd) %>% 
      rename( R = R_overall,
              `R 2.5% CI` = `2.5%`,
              `R 97.5% CI` = `97.5%`),
    by = c("Trait", "Trait_Sex")) %>% 
  left_join(
    R_diff_transcript_summarised %>%
      select(-sd) %>% 
      rename(`Delta R` = R_diff,
             `Delta R 2.5% CI` = `2.5%`,
             `Delta R 97.5% CI` = `97.5%`),
    by = c("Trait", "Trait_Sex")) %>% 
  
  left_join(evidence_ratios_transcriptome %>% 
              mutate(Trait_Sex = case_when(Trait_Sex == "Traits measured in females" ~ "Female",
                                           Trait_Sex == "Traits measured in males" ~ "Male")) %>% 
              select(-c(entrez_id, chromosome, chromosome_numeric, position)) %>% 
              rename(`R ER` = Overall_Evid_Ratio, `R female ER` = Female_Evid_Ratio,
                     `R male ER` = Male_Evid_Ratio, `Delta R ER` = Diff_Evid_Ratio,
                     `Sexual concordance ER` = Concord_Evid_Ratio)) %>% 
  rename(`Sex trait was measured in` = Trait_Sex) %>% 
    select(Trait, gene_name, gene_symbol, `Sex trait was measured in`, R, `R 2.5% CI`, `R 97.5% CI`, `R ER`, 
           `R female`, `R female 2.5% CI`, `R female 97.5% CI`, `R female ER`,
           `R male`, `R male 2.5% CI`, `R male 97.5% CI`, `R male ER`,
           `Delta R`, `Delta R 2.5% CI`, `Delta R 97.5% CI`, `Delta R ER`, `Sexual concordance ER`) %>% 
  mutate(across(c(5,6,7,9,10,11,13,14,15,17,18,19), ~round(.x, 3)),
         across(ends_with("ER"), ~ round(.x, 4))) 

write_csv(all_transcriptome_summary, "data/transcriptome_output/Supplementary_dataset_2.csv")

my_data_table(all_transcriptome_summary)

\(~\)

Estimate \(\overline{|R|}\) in each sex

\(~\)

The model

Fit the model to test whether \(\overline{|R|}\) depends on the sex fitness and trait values were measured in:

# fit the model

median_R_transcriptome_model <- 
  brm(absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1|Trait),
      family = brmsfamily(family = "Gamma"), # gamma is appropriate for the half-normal distribution created by taking the absolute
      data = R_summarised_transcriptome,
      prior = c(prior(normal(-2.2, 1), class = Intercept),
                prior(exponential(1), class = sd),
                prior(exponential(1), class = shape),
                prior(normal(0, 0.25), class = b)),
      warmup = 2000, iter = 6000,
      seed = 1, cores = 4, chains = 4,
      control = list(adapt_delta = 0.8, max_treedepth = 10),
      file = "fits/median_R_transcriptome_model")

print(median_R_transcriptome_model)

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1 | Trait) 
   Data: R_summarised_transcriptome (Number of observations: 57072) 
  Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
         total post-warmup draws = 16000

Group-Level Effects: 
~Trait (Number of levels: 14268) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.41      0.00     0.41     0.42 1.00     2388     5240

Population-Level Effects: 
                              Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                        -2.72      0.00    -2.72    -2.71 1.00
Fitness_SexMale                  -0.01      0.00    -0.01    -0.00 1.00
Trait_SexMale                    -0.11      0.00    -0.11    -0.10 1.00
Fitness_SexMale:Trait_SexMale     0.03      0.00     0.03     0.04 1.00
                              Bulk_ESS Tail_ESS
Intercept                         1945     4332
Fitness_SexMale                  15332    12393
Trait_SexMale                    15791    12963
Fitness_SexMale:Trait_SexMale    15058    12576

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     1.65      0.00     1.65     1.66 1.00    23558    11868

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

invisible(gc())

To check how does the gamma distribution fares at capturing the shape of the data, we use the model to recapitulate the data it was trained on. Given the shape of the predicted data is similar to the real data, there’s a high chance we’ve used an appropriate distribution family.

pp_check(median_R_transcriptome_model)

Past versions of unnamed-chunk-37-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11

\(~\)

Get model predictions. Once again, we only do this once and save for later convenience.

new_data <- expand_grid(Fitness_Sex = c("Female", "Male"),
                        Trait_Sex = c("Female", "Male"))

if(!file.exists("data/transcriptome_output/R_transcriptome_fitted.csv")){
R_transcriptome_fitted <-
  fitted(median_R_transcriptome_model, newdata = new_data, re_formula = NA, summary = F) %>% 
  as.data.frame() %>% 
  rename(FemaleFitness_FemaleTrait = V1, FemaleFitness_MaleTrait = V2, 
         MaleFitness_FemaleTrait = V3, MaleFitness_MaleTrait = V4) %>% 
  as_tibble() %>% 
  mutate(percent_diff_female = 
           ((MaleFitness_FemaleTrait - FemaleFitness_FemaleTrait) / FemaleFitness_FemaleTrait)*100,
         percent_diff_male = 
           ((MaleFitness_MaleTrait - FemaleFitness_MaleTrait)/ FemaleFitness_MaleTrait)*100) %>% 
  pivot_longer(cols = everything(), names_to = "Parameter", values_to = "R_mean") %>% 
  mutate(Fitness_Sex = case_when(str_detect(Parameter, "FemaleFitness") ~ "Female",
                                 str_detect(Parameter, "MaleFitness") ~ "Male"),
         Trait_Sex = case_when(str_detect(Parameter, "FemaleTrait") ~ "Female",
                               str_detect(Parameter, "MaleTrait") ~ "Male",
                               str_detect(Parameter, "percent_diff_female") ~ "Female",
                               str_detect(Parameter, "percent_diff_male") ~ "Male"))

write_csv(R_transcriptome_fitted, "data/transcriptome_output/R_transcriptome_fitted.csv")
} else R_transcriptome_fitted <- read_csv("data/transcriptome_output/R_transcriptome_fitted.csv")

\(~\)

Complete Figure 2

Create panels b and d

R_t1 <-
  R_transcriptome_fitted %>% 
  filter(!str_detect(Parameter, "percent")) %>% 
  ggplot(aes(x = R_mean, y = Trait_Sex, fill = Fitness_Sex)) + #+ y = Trait_Sex)) + #fill = Fitness_Sex)) +
  stat_slab(alpha = 0.9, shape = 21) +#,
  labs(x = expression("Mean |"* italic(R) * "| for gene expression traits"),
       y = "Phenotyped sex",
       fill = "Sex under selection") +
  scale_fill_manual(values = c(carto_pal(7, "Purp")[5], carto_pal(7, "Peach")[5])) +
  coord_cartesian(xlim = c(0.045, 0.10)) +
  theme_minimal() + 
  theme(panel.background = element_rect(fill='transparent'),
        #panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA),
        legend.position = "bottom",
        text = element_text(size=13))

R_t2 <-
  R_transcriptome_fitted %>% 
  filter(str_detect(Parameter, "percent")) %>%
  ggplot(aes(x = R_mean, y = Trait_Sex)) + #, #y = Trait_Sex, fill = Trait_Sex)) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, fill = carto_pal(7, "Emrld")[2],
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.2) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals + # width indicates the uncertainty intervals: here we have 66% and 95% intervals+
  #scale_fill_manual(values = c(carto_pal(7, "Peach")[5], carto_pal(7, "Peach")[5], carto_pal(7, "Purp")[5])) +
  geom_vline(xintercept = 0, linetype = 2, linewidth = 1.2) +
  coord_cartesian(xlim = c(-5, 45)) +
  xlab(expression("% sex difference in |"*italic(R)* "|")) +
  ylab("Phenotyped sex") +
  theme_minimal() + 
  theme(panel.background = element_rect(fill='transparent'),
        #panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA),
        legend.position = "none",
        text = element_text(size=13))

Produce the complete figure

left_col <- (R_p1 / R_t1)

right_col <- (R_p2 / R_t2)

wrap_plots(left_col, right_col) & 
  plot_annotation(tag_levels = "a") 

Past versions of unnamed-chunk-40-1.png

Version	Author	Date
ad196cf	tomkeaney	2024-04-11

Figure 2. a the posterior distribution for |\(\overline{R}\)| across organismal level traits, the mean expected response to selection (absolute value) in females and males and split by the sex the trait was measured in. b shows |\(\overline{R}\)| across gene-expression traits. |\(\overline{R}\)| is expressed in trait standard deviations. c and d show the percentage increase or decrease in |\(\overline{R}\)| through males compared to through females for the organismal level and gene expression traits, respectively. Points indicate the estimated mean with associated 66 and 95% credible intervals (where wide enough to be visible).

Calculating stats for the results section

R_transcriptome_fitted %>% 
  group_by(Parameter) %>% 
  median_qi(R_mean)

# A tibble: 6 × 7
  Parameter                  R_mean  .lower  .upper .width .point .interval
  <chr>                       <dbl>   <dbl>   <dbl>  <dbl> <chr>  <chr>    
1 FemaleFitness_FemaleTrait  0.0661  0.0656  0.0666   0.95 median qi       
2 FemaleFitness_MaleTrait    0.0593  0.0588  0.0597   0.95 median qi       
3 MaleFitness_FemaleTrait    0.0655  0.0650  0.0660   0.95 median qi       
4 MaleFitness_MaleTrait      0.0608  0.0603  0.0612   0.95 median qi       
5 percent_diff_female       -0.890  -1.47   -0.316    0.95 median qi       
6 percent_diff_male          2.55    1.96    3.13     0.95 median qi       

Session information

sessionInfo()

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

Matrix products: default


locale:
[1] LC_COLLATE=English_Germany.utf8  LC_CTYPE=English_Germany.utf8   
[3] LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C                    
[5] LC_TIME=English_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bigsnpr_1.12.4    bigstatsr_1.5.12  ggnewscale_0.4.9  ggrepel_0.9.4    
 [5] ggtext_0.1.2      broom_1.0.5       bayestestR_0.13.1 tidybayes_3.0.6  
 [9] brms_2.20.4       Rcpp_1.0.11       patchwork_1.1.3   DT_0.30          
[13] kableExtra_1.3.4  rcartocolor_2.1.1 MetBrewer_0.2.0   lubridate_1.9.3  
[17] forcats_1.0.0     stringr_1.5.0     dplyr_1.1.3       purrr_1.0.2      
[21] readr_2.1.4       tidyr_1.3.0       tibble_3.2.1      ggplot2_3.4.4    
[25] tidyverse_2.0.0   workflowr_1.7.1  

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2       rstudioapi_0.15.0    jsonlite_1.8.7      
  [4] magrittr_2.0.3       farver_2.1.1         rmarkdown_2.25      
  [7] fs_1.6.3             vctrs_0.6.4          base64enc_0.1-3     
 [10] webshot_0.5.5        htmltools_0.5.6.1    distributional_0.3.2
 [13] curl_5.1.0           sass_0.4.7           StanHeaders_2.26.28 
 [16] bslib_0.5.1          htmlwidgets_1.6.2    plyr_1.8.9          
 [19] zoo_1.8-12           cachem_1.0.8         commonmark_1.9.0    
 [22] whisker_0.4.1        igraph_1.5.1         iterators_1.0.14    
 [25] mime_0.12            lifecycle_1.0.4      pkgconfig_2.0.3     
 [28] colourpicker_1.3.0   Matrix_1.6-1.1       R6_2.5.1            
 [31] fastmap_1.1.1        shiny_1.7.5.1        digest_0.6.33       
 [34] colorspace_2.1-0     ps_1.7.5             rprojroot_2.0.3     
 [37] crosstalk_1.2.0      labeling_0.4.3       fansi_1.0.5         
 [40] timechange_0.2.0     httr_1.4.7           abind_1.4-5         
 [43] compiler_4.3.1       rngtools_1.5.2       bit64_4.0.5         
 [46] doParallel_1.0.17    withr_2.5.1          backports_1.4.1     
 [49] inline_0.3.19        shinystan_2.6.0      bigassertr_0.1.6    
 [52] QuickJSR_1.0.7       pkgbuild_1.4.2       gtools_3.9.4        
 [55] loo_2.6.0            tools_4.3.1          httpuv_1.6.11       
 [58] threejs_0.3.3        glue_1.6.2           callr_3.7.3         
 [61] nlme_3.1-162         promises_1.2.1       cmdstanr_0.6.0      
 [64] bigparallelr_0.3.2   gridtext_0.1.5       grid_4.3.1          
 [67] checkmate_2.2.0      getPass_0.2-2        reshape2_1.4.4      
 [70] generics_0.1.3       gtable_0.3.4         bigsparser_0.6.1    
 [73] tzdb_0.4.0           flock_0.7            data.table_1.14.8   
 [76] hms_1.1.3            xml2_1.3.5           utf8_1.2.3          
 [79] foreach_1.5.2        pillar_1.9.0         ggdist_3.3.0        
 [82] markdown_1.10        vroom_1.6.4          posterior_1.4.1     
 [85] later_1.3.1          lattice_0.21-8       bit_4.0.5           
 [88] tidyselect_1.2.0     miniUI_0.1.1.1       knitr_1.44          
 [91] git2r_0.32.0         arrayhelpers_1.1-0   gridExtra_2.3       
 [94] V8_4.4.0             svglite_2.1.2        stats4_4.3.1        
 [97] xfun_0.40            bridgesampling_1.1-2 matrixStats_1.0.0   
[100] rstan_2.32.3         stringi_1.7.12       yaml_2.3.7          
[103] evaluate_0.22        codetools_0.2-19     cli_3.6.1           
[106] RcppParallel_5.1.7   shinythemes_1.2.0    xtable_1.8-4        
[109] systemfonts_1.0.5    munsell_0.5.0        processx_3.8.2      
[112] jquerylib_0.1.4      coda_0.19-4          svUnit_1.0.6        
[115] parallel_4.3.1       rstantools_2.3.1.1   ellipsis_0.3.2      
[118] prettyunits_1.2.0    doRNG_1.8.6          dygraphs_1.1.1.6    
[121] bayesplot_1.10.0     Brobdingnag_1.2-9    viridisLite_0.4.2   
[124] mvtnorm_1.2-3        scales_1.3.0         xts_0.13.1          
[127] insight_0.19.6       crayon_1.5.2         rlang_1.1.1         
[130] cowplot_1.1.1        rvest_1.0.3          shinyjs_2.1.0